Discussion of “Geodesic Monte Carlo on Embedded Manifolds” Comment: Connections and Extensions

We welcome this paper of Byrne and Girolami [BG]; it breathes even more life into the emerging area of hybrid Monte Carlo Markov chains by introducing original tools for dealing with Monte Carlo simulations on constrained spaces such as manifolds. We begin our comment with a bit of history. Using geodesics to sample from the uniform distribution on Stiefel manifold was proposed by Asimov (1985) in his work on the Grand Tour for exploratory data analysis. For data x1, x2,..., xn in Rp, it is natural to inspect low dimensional projections γx1, γx2,..., γxn for γ: Rp − → Rk. In the [BG] paper the authors have a space of k-frames in Rp, called Vk,p. If one chooses γ at random from this space, the views would be too ‘disconnected ’ or ‘jerky ’ for human observers. A better tactic turned out to be to choose a few γi, 1 ≤ i ≤ L, at random and then moving smoothly from γi to γj by available closed form geodesics. While in a historical mode, we point to the little known papers of McLachlan and Quispel (2003) and more recent papers by Betancourt on hybrid Monte Carlo (Betancourt, 2013). Discrete Hamiltonian Dynamics The paper of [BG] uses Hamiltonian dynamics to move around on a manifold in an intelligent way to get proposals for the Metropolis algorithm. There are also many problems where samples are needed for constrained discrete spaces. These include sampling contingency tables with given row and column sums as in Diaconis and Sturmfels (1998). We recently encountered the following prob-lem in a quantum physics context (Chatterjee and Diaconis, 2013). Consider boxes labeled (1, 2, 3,...). Drop N balls into these boxes according to Bose-Einstein allocation resulting in Ni balls in box labeled i. Interest is on samples conditional on iNii 2 = E. This is a discrete version of the author’s sam-pling from simplices and spheres. We do not currently have discrete versions of Hamiltonian dynamics apart from numerical schemes (leapfrog) that are used to solve the resulting differential equations as proposed by Neal (2011). In contrast, [BG] compute the dynamics by splitting up the Hamiltonian into two analyti-cally solvable parts. We wonder whether the author’s can suggest adaptions of their ideas to the discrete framework